Script introduction
Find in the following lines the code corresponding to the analysis of the article: Heterogeneity in the isolation of patches may be essential for the action of metacommunity mechanisms.
All the suplementary material images, original scripts and used functions can be found in:
Article DOI:
Direct link:
To properly follow the following steps we recommend to read the entire work to better contextualise the use and limitations of the current codes and functions.
The scripts can be found in the following repository
Link to the function: https://github.com/ Link to supplementary material:
# Load basic functions to run all the simulations
library(dplyr); library(tibble); library(ggplot2); library(gridExtra)Main developedfunctoins
To develpe and carry all the analysis that readers can found in the previously cited paper we developed two main functions that define our two main analysis:
1. Pondscape diversity along dispersal and centrality gradients (mainly providing Figure 3A and B columns).
2. Pondscape randomization and diveristy comparisons (mainly providing Figure 3C and D colums)
For each one of these analysis we provide in the repository the scripts developing the corresponding functions: 1) FEE_IterationModel_Ponderful.R and 2) FEE_Random_landscape_Ponderful.R.
In these scripts there is a detailed explanations of what each line is conducting and which is the specific objective deboted to each function. We begun the analysis by charging these two functions as well as the .RData archive where all the necessary information is contained (xy coordinates of each water body from each one of the studied pondscapes).
# Charging functoins
source("FEE_Random_landscape_Ponderful.R")
source("FEE_IterationModel_Ponderful.R")
# We load all the necessary information for each studied pondscape
load("Espacio_Retromed.RData")
CHARCOS_ROCHA <- read.delim("CHARCOS_ROCHA.txt")
Rocha<-CHARCOS_ROCHA
Rocha <- cbind(Rocha,0,1,100)Lines of analysis for each Pondscape
For each podnscape we run the same set of code lines where we define the set of requirements fo the two main functions utilisation. The following lines can also be found in the script named Lines_Analysis.R contained in the main github repository of this project (see above).
##Albera
# Seting elements for each function
Albera<-cbind(ALBERA,100)
Albera<-Albera[,-4]
colnames(Albera)[8]<-"J"
Albera<-Albera[,c(1:3, 5,7,6,4,8)]
M.distance.Albera<-as.matrix(dist(Albera[,2:3]))
Meta.comm.Albera<-Meta.comm.Albera
J.Albera<-Albera[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Albera$grado
Cc<-Albera$cc
Bt<-Albera$bc
N<-nrow(M.distance.Albera)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Albera<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Albera, m.pool=0.01, Js=J.Albera, id.module=Albera[,7],
M.dist=M.distance.Albera, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Albera,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Albera))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Albera<-rbind(out.d50.obs.landscape.Albera, out.t2)
}
out.d50.obs.landscape.Albera
save.image("~/Dropbox/H2020/RETROMED/Frontiers_Diciembre_2022/Frontiers_Respuestas/Retromed_H2020_febrero_2023.RData")
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Albera_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Albera, m.pool=0.01,
Js=Albera[,8],
id.module=Albera[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Albera[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Albera,Lottery=F, it=0,
prop.dead.by.it=0,id.obs=1:nrow(Albera))##Giara
# Seting elements for each function
Giara<-cbind(GIARA,100)
Giara<-Giara[,-4]
colnames(Giara)[8]<-"J"
Giara<-Giara[,c(1:3, 5,7,6,4,8)]
M.distance.Giara<-as.matrix(dist(Giara[,2:3]))
Meta.comm.Giara<-Meta.comm
J.Giara<-Giara[,8]
#D50s<-seq(50,3000,,50)
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Giara$grado
Cc<-Giara$cc
Bt<-Giara$bc
N<-nrow(M.distance.Giara)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Giara<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Giara, m.pool=0.01, Js=J.Giara, id.module=Giara[,7],
M.dist=M.distance.Giara, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Giara,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Giara))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Giara<-rbind(out.d50.obs.landscape.Giara, out.t2)
}
out.d50.obs.landscape.Giara
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Giara_febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Giara, m.pool=0.01, Js=Giara[,8],
id.module=Giara[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Giara[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Giara,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Giara))##Guils
# Seting elements for each function
Guils<-cbind(GUILS,100)
Guils<-Guils[,-4]
colnames(Guils)[8]<-"J"
Guils<-Guils[,c(1:3, 5,7,6,4,8)]
M.distance.Guils<-as.matrix(dist(Guils[,2:3]))
Meta.comm.Guils<-Meta.comm.Albera
J.Guils<-Guils[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Guils$grado
Cc<-Guils$cc
Bt<-Guils$bc
N<-nrow(M.distance.Guils)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Guils<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Guils, m.pool=0.01, Js=J.Guils, id.module=Guils[,7],
M.dist=M.distance.Guils, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Guils,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Guils))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Guils<-rbind(out.d50.obs.landscape.Guils, out.t2)
}
out.d50.obs.landscape.Guils
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Guils_febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 1000,Meta.pool=Meta.comm.Guils, m.pool=0.01,
Js=Guils[,8],
id.module=Guils[,7],
D50.min=1, D50.max=2000,m.max=1,M.X.Y=Guils[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Guils,
Lottery=F, it=0, prop.dead.by.it=0,
id.obs=1:nrow(Guils))##PNAE
# Seting elements for each function
Pnae<-cbind(PNAE,100)
Pnae<-Pnae[,-4]
colnames(Pnae)[8]<-"J"
Pnae<-Pnae[,c(1:3, 5,7,6,4,8)]
M.distance.Pnae<-as.matrix(dist(Pnae[,2:3]))
Meta.comm.Pnae<-Meta.comm
J.Pnae<-Pnae[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Pnae$grado
Cc<-Pnae$cc
Bt<-Pnae$bc
N<-nrow(M.distance.Pnae)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Pnae<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Pnae, m.pool=0.01, Js=J.Pnae, id.module=Pnae[,7],
M.dist=M.distance.Pnae, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Pnae,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Pnae))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Pnae<-rbind(out.d50.obs.landscape.Pnae, out.t2)
}
out.d50.obs.landscape.Pnae
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Pnae_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Pnae, m.pool=0.01, Js=Pnae[,8],
id.module=Pnae[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Pnae[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Pnae,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Pnae))##Vila nova
# Seting elements for each function
Vilavona<-cbind(VILAVONA,100)
Vilavona<-Vilavona[,-4]
colnames(Vilavona)[8]<-"J"
Vilavona<-Vilavona[,c(1:3, 5,7,6,4,8)]
M.distance.Vilavona<-as.matrix(dist(Vilavona[,2:3]))
Meta.comm.Vilavona<-Meta.comm.Albera
J.Vilavona<-Vilavona[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Vilavona$grado
Cc<-Vilavona$cc
Bt<-Vilavona$bc
N<-nrow(M.distance.Vilavona)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Vilavona<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Vilavona, m.pool=0.01, Js=J.Vilavona, id.module=Vilavona[,7],
M.dist=M.distance.Vilavona, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Vilavona,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Vilavona))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Vilavona<-rbind(out.d50.obs.landscape.Vilavona, out.t2)
}
out.d50.obs.landscape.Vilavona
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Vilavona_Febrero_2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Vilavona,
m.pool=0.01, Js=Vilavona[,8],
id.module=Vilavona[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Vilavona[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1,
comm.fixed=Meta.comm.Vilavona,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Vilavona))##Rocha
# Seting elements for each function
colnames(Rocha)[6:8]<-c("bc","modules", "J")
Rocha[1:10,7]<-2
M.distance.Rocha<-as.matrix(dist(Rocha[,2:3]))
Meta.comm.Rocha<-Meta.comm.Albera
J.Rocha<-Rocha[,8]
D50s<-10^seq(log10(10),log10(3000),,50)
Gr<-Rocha$grado
Cc<-Rocha$cc
Bt<-Rocha$bc
N<-nrow(M.distance.Rocha)
# Utilisation of the "FEE_IterationModel_Ponderful.R" functions
out.d50.obs.landscape.Rocha<-NULL
for (d in D50s){
out.t<-iteration.model(replicas=112, Meta.pool=Meta.comm.Rocha, m.pool=0.01, Js=J.Rocha, id.module=Rocha[,7],
M.dist=M.distance.Rocha, D50=d, m.max=1,
id.fixed=NULL, D50.fixed=1000, m.max.fixed=1, comm.fixed=Meta.comm.Rocha,
Lottery=F, it=2000, prop.dead.by.it=0.05, id.obs=1:ncol(M.distance.Rocha))
out.t<-resume.out(out.t)
S<-out.t[[1]][11:(10+N)]
Be.all<-out.t[[1]][(10+N+1):(10+N+N)]
gamma<-out.t[[1]][10]
print(c(length(d),length(S), length(Be.all), length(gamma),length(Bt), length(Gr), length(Cc)))
out.t2<-cbind(d,S, Be.all, gamma,Bt, Gr, Cc)
print(tail(out.t2))
out.d50.obs.landscape.Rocha<-rbind(out.d50.obs.landscape.Rocha, out.t2)
}
out.d50.obs.landscape.Rocha
# Utilisation of the "FEE_Random_landscape_Ponderful.R" functions
# sequence of D50
random_vs_real_D50.Rocha_febrero2023<-H2020_Random_Landscape_local_ponds.D50(n.random = 140,Meta.pool=Meta.comm.Rocha, m.pool=0.01, Js=Rocha[,8],
id.module=Rocha[,7],
D50.min=1, D50.max=2000, m.max=1,M.X.Y=Rocha[,2:3],
id.fixed=0, D50.fixed=100, m.max.fixed=1, comm.fixed=Meta.comm.Rocha,
Lottery=F, it=0, prop.dead.by.it=0, id.obs=1:nrow(Rocha))Lines of ploting and results representation for each Pondscape
Once all the results have been produced and all the simulations carried we extract the specific outputs of alpha diversity and betadiversity in order to plot them accrodingly and to represent the article figures. The following lines can also be found in the script named Lines_PLots.R contained in the main github repository of this project (see above).
Figure 1
xy.ALB <- read.table(file = "Pondscapes_loc/CHARCOS_ALBERA.txt",sep = "\t", dec =",",header = T)
xy.GIA <- read.table(file = "Pondscapes_loc/CHARCOS_GIARA.txt",sep = "\t", dec =",",header = T)
xy.GUI <- read.table(file = "Pondscapes_loc/CHARCOS_GUILS.txt",sep = "\t", dec =",",header = T)
xy.PNAE <- read.table(file = "Pondscapes_loc/CHARCOS_PNAE.txt",sep = "\t", dec =",",header = T)
xy.VMIL <- read.table(file = "Pondscapes_loc/CHARCOS_VILANOVAMILFONTES.txt",sep = "\t", dec =",",header = T)
xy.ROC <- read.table(file = "Pondscapes_loc/CHARCOS_ROCHA.txt",sep = "\t", dec =".",header = T)
# Check structure and correct PNAE
str(xy.ALB)
str(xy.GIA)
str(xy.GUI)# Remember that they are repeated
xy.GUI<- xy.GUI[-which(duplicated(xy.GUI[,2:3])==T),]
str(xy.PNAE)# Remember that they are repeated
xy.PNAE<- xy.PNAE[-which(duplicated(xy.PNAE[,2:3])==T),]
str(xy.VMIL)
str(xy.ROC)
c(nrow(xy.ALB),nrow(xy.GIA),nrow(xy.GUI),nrow(xy.PNAE),nrow(xy.VMIL),nrow(xy.ROC))
c(542,160,86,255,146)
# Creation of minimum spanning tree (one link per pond)
GG_ALB <- mst(as.matrix(dist(xy.ALB)))
GG_GIA <- mst(as.matrix(dist(xy.GIA)))
GG_GUI <- mst(as.matrix(dist(xy.GUI)))
GG_PNAE <- mst(as.matrix(dist(xy.PNAE)))
GG_VMIL <- mst(as.matrix(dist(xy.VMIL)))
GG_ROC <- mst(as.matrix(dist(cbind(xy.ROC[,2],xy.ROC[,3]))))
# List creation
xy.values <- list(xy.ALB[,2:3],xy.PNAE[,2:3],xy.GIA[,2:3],xy.VMIL[,1:2], xy.GUI[,2:3], xy.ROC[,2:3])
GG.values <- list(GG_ALB,GG_PNAE, GG_GIA, GG_VMIL, GG_GUI, GG_ROC)
Pondscapes_names <- c("Albera","PNAE","Giara", "Vila Nova","Guils", "Rocha")
# Nice colours? CUNILLERA_palette is what you need
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R") # Also found in https://github.com/Cunillera-Montcusi/Cunillera_Pallete
library(sna)
gg_plot_network <- list()
for (u in 1:length(xy.values)) {
detach("package:sna", unload = TRUE)
library(igraph)
# In order to plot the network in a ggplot way do the follwing
weight_degree <- Pondscapes_outputs[[u]][which(Pondscapes_outputs[[u]][,1]==max_diff[u]),7]
detach("package:igraph", unload = TRUE)
library("sna")
n<- network(GG.values[[u]], directed=F, diag=F)
n %v% "family" <- weight_degree # Family is an standard name...
# ALL THESE LINES MUST BE RUNNED!!!
#from here_____________________________________________________________________
gg_plot_network[[u]] <- ggplot(n, layout=as.matrix(xy.values[[u]]),
aes(x = x, y = y, xend = xend, yend = yend))+
geom_edges(color = "grey30", size=0.5) +
geom_nodes(aes(fill=family, size=family) ,color ="black" ,shape=21, alpha=.55)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(x="",y="",title = Pondscapes_names[[u]])+
theme(axis.line = element_blank(),
plot.title = element_text(size=35, colour = "black",hjust = 0.5),
axis.text = element_blank(),
axis.ticks = element_blank(),
legend.position = "none",
panel.background = element_blank())
legend_plots<- get_legend(ggplot(n, layout=as.matrix(xy.values[[u]]),
aes(x = x, y = y, xend = xend, yend = yend))+
geom_edges(color = "black", size=0.5) +
geom_nodes(aes(fill=family) ,color ="grey50" ,shape=21, alpha=.55)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Weighted degree",
guide = guide_colorbar(ticks = T,
barheight = 5,
barwidth = 1))+
theme_classic()+
theme(legend.direction = "vertical",
legend.box="vertical"))
}
Figure 3 settings
# Nice colours? CUNILLERA_palette is what you need
source("C:/Users/David CM/Dropbox/DAVID DOC/LLAM al DIA/CUNILLERA_palette.R")
Pondscapes_outputs <- list(out.d50.obs.landscape.Albera,
out.d50.obs.landscape.Pnae,
out.d50.obs.landscape.Giara,
out.d50.obs.landscape.Vilavona,
out.d50.obs.landscape.Guils,
out.d50.obs.landscape.Rocha)
Pondscapes__rand_outputs <- list(random_vs_real_D50.Albera_Febrero_2023,
random_vs_real_D50.Pnae_Febrero_2023,
random_vs_real_D50.Giara_febrero_2023,
random_vs_real_D50.Vilavona_Febrero_2023,
random_vs_real_D50.Guils_Febrero_2023,
random_vs_real_D50.Rocha_febrero2023)
Pondscapes_names <- c("Albera", "PNAE", "Giara", "Vila Nova", "Guils","Rocha")
final_plot_v1 <- list()
final_plot_v2 <- list()
supp_CV_plot <- list()Figure 3 Albera (first row)
#Albera________________________________________________________________________________________________________####
p=1
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=30, q=3, d50=150))
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=32, q=3, d50=100))
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>% ggplot(aes(d,S))+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(title="A)",subtitle = Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(40), xmax = log10(250), ymin = 20, ymax = 37)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=500))
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-3, d50=150))
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>% arrange(desc(d)) %>% ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(title="B)",subtitle = " ", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(40), xmax = log10(250), ymin = .5, ymax = .8)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = "C)", subtitle = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title="D)", subtitle = " ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- data.frame()
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Albera",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 PNAE (second row)
#PNAE________________________________________________________________________________________________________####
p=2
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.01))
ii.max.g<-which(a[,6]>quantile(a[,6],.99))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=0, Smax=30, q=3, d50=500))
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=0, Smax=20, q=3, d50=100))
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(520), xmax = log10(2900), ymin = 1, ymax = 22)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-1, d50=200) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(510), xmax = log10(3100), ymin = .7, ymax = .99)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+#scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+#scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "PNAE",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 Giara (third row)
#Giara________________________________________________________________________________________________________####
p=3
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=0, Smax=25, q=3, d50=500) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=35, q=1, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle = Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(900), xmax = log10(3200), ymin = 0, ymax = 22)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=500) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=10, Bmax=1, q=-3, d50=50) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(40), xmax = log10(250), ymin = .42, ymax = .74)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+#scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+#scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Giara",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 Vila nova (foruth row)
#Vilanova________________________________________________________________________________________________________####
p=4
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=30, q=3, d50=200) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=32, q=3, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(400), xmax = log10(3000), ymin = 1, ymax = 27)
#### BETA
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-3, d50=150) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for Beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(500), xmax = log10(3000), ymin = .53, ymax = .99)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Vila Nova",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 Guils (fifth row)
#Guils________________________________________________________________________________________________________####
p=5
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=25, q=3, d50=500) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=35, q=1, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(600), xmax = log10(3000), ymin = 0.5, ymax = 22.5)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-1, d50=200) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(650), xmax = log10(3100), ymin = .55, ymax = .99)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title=" ")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Guils",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 Rocha (sixth row)
#ROCHA ________________________________________________________________________________________________________####
p=6
a<- Pondscapes_outputs[[p]]
############# Figures Albera
#########################################3
### Plot results: alpha and beta diversity in species` dispersal gradient
### also considering the gradient in local ponds isolation
# Calculation of the minumum and maximum quantiles
ii.min.g<-which(a[,6]<quantile(a[,6],.05))
ii.max.g<-which(a[,6]>quantile(a[,6],.95))
# Minimum quantile line
d<-a[ii.min.g,1]
S<-a[ii.min.g,2]
gg<-a[ii.min.g,6]
cc<-a[ii.min.g,7]
B<-a[ii.min.g,3]
# Maximum quantile line
d2<-a[ii.max.g,1]
S2<-a[ii.max.g,2]
gg2<-a[ii.max.g,6]
cc2<-a[ii.max.g,7]
B2<-a[ii.max.g,3]
### ALPHA - S
#Function to descrive the minimum curve
M.hill<-nls(S~S0+(Smax-S0)*(d^q)/(d50^q+d^q), start = list(S0=4, Smax=30, q=3, d50=200) )
summary(M.hill)
ph<-coefficients(M.hill)
f.hill<-function(x) ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3])
#Function to descrive the maximum curve
M.hill.2<-nls(S2~S0+(Smax-S0)*(d2^q)/(d50^q+d2^q), start = list(S0=4, Smax=32, q=3, d50=100) )
summary(M.hill.2)
ph2<-coefficients(M.hill.2)
f.hill.2<-function(x) ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3])
# Differential delta between the two lines and small plot of it
f.hill.delta<-function(x)(ph2[1]+(ph2[2]-ph2[1])*(x^ph2[3])/(ph2[4]^ph2[3]+x^ ph2[3]))/(ph[1]+(ph[2]-ph[1])*(x^ph[3])/(ph[4]^ph[3]+x^ ph[3]))
S.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Alpha ratio",subtitle = "Alpha ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for alpha
S.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,S)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.2 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle=Pondscapes_names[p], y="Alpha diversity", x=NULL)+
geom_function(fun=f.hill, col="white", size=1.8)+
geom_function(fun=f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=f.hill.2, col="white", size=1.8)+
geom_function(fun=f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
S.o.Albera<-S.o.Albera+annotation_custom(ggplotGrob(S.ratio), xmin = log10(300), xmax = log10(2000), ymin = 2, ymax = 19)
#### Beta
#Function to descrive the minimum curve
B.M.hill<-nls(B~B0+(Bmax-B0)*(d^q)/(d50^q+d^q), start = list(B0=4.396e-01, Bmax=9.179e-01, q=-2.56, d50=250) )
summary(B.M.hill)
B.ph<-coefficients(B.M.hill)
B.f.hill<-function(x) B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3])
#Function to descrive the maximum curve
B.M.hill.2<-nls(B2~B0+(Bmax-B0)*(d2^q)/(d50^q+d2^q), start = list(B0=0, Bmax=10, q=-3, d50=150) )
summary(B.M.hill.2)
B.ph2<-coefficients(B.M.hill.2)
B.f.hill.2<-function(x) B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3])
# Differential delta between the two lines and small plot of it
B.f.hill.delta<-function(x)(B.ph2[1]+(B.ph2[2]-B.ph2[1])*(x^B.ph2[3])/(B.ph2[4]^B.ph2[3]+x^B.ph2[3]))/(B.ph[1]+(B.ph[2]-B.ph[1])*(x^B.ph[3])/(B.ph[4]^B.ph[3]+x^B.ph[3]))
B.ratio<-ggplot(data.frame(a))+aes(d)+geom_blank()+
geom_function(fun=B.f.hill.delta, col="red")+
theme_classic()+labs(x="dispersal",y="Beta ratio",subtitle = "Beta ratio")+
scale_x_log10()+
theme(axis.text=element_text(size=5),axis.title=element_text(size=6,face="bold"))
# Plot for beta
B.o.Albera<-data.frame(a) %>%arrange(desc(d)) %>%ggplot()+aes(d,Be.all)+
#geom_point(aes(fill=Gr),shape=21, alpha=0.2,size=3, colour="black")+
geom_jitter(aes(fill=Gr,colour=Gr),height =0.02 ,shape=21, alpha=0.2,size=2)+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
scale_color_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F)+
labs(subtitle="", y="Beta diversity", x=NULL)+
geom_function(fun=B.f.hill, col="white", size=1.8)+
geom_function(fun=B.f.hill, col="black", linetype="dashed", size=1.5)+
geom_function(fun=B.f.hill.2, col="white", size=1.8)+
geom_function(fun=B.f.hill.2, col="black", linetype="solid", size=1.5)+
xlab("Dispersal (d50)")+
scale_x_log10()+theme_bw()+theme(legend.position = "none")
B.o.Albera<-B.o.Albera+annotation_custom(ggplotGrob(B.ratio), xmin = log10(500), xmax = log10(2500), ymin = .5, ymax = .95)
#Legend extraction of the gradient of Closeness for both alpha and beta (it is the same for both)
legend <- get_legend(ggplot(data.frame(a))+aes(d,Be.all, fill=Gr)+
geom_point(shape=21, alpha=0.1,size=3, colour="black")+
scale_fill_CUNILLERA(palette = "LGTBI", discrete = F, reverse = F, name="Degree")+
scale_x_log10()+theme_bw())
# Obtention of Z-plots responding to the differential between alpha diversity and beta diversity
#Check function "plotea.null.landscape.D50" for more information
b<-data.frame(Pondscapes__rand_outputs[[p]])
b<-mutate(b, Z.S=(mean.S.real.-mean.S.null.)/sd.S.null.)
b$Z.S[which(b$Z.S==Inf)] <- 0
b<-mutate(b, Z.Bet=(mean.Be.all.real.-mean.Be.all.null.)/sd.Be.all.null.)
b<-mutate(b, CV.oe.S=(sd.S.real./mean.S.real.-sd.S.null./mean.S.null.))
b<-mutate(b, CV.oe.Bet=(sd.Be.all.real./mean.Be.all.real.-sd.Be.all.null./mean.Be.all.null.))
b<-mutate(b, Z.Gama=(Gamma.real-Gamma.null)/sd(Gamma.null))
S.oe<-ggplot(b, aes(D50,Z.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
labs(y="Z-Alpha", x="Dispersal (d50)", size=0.25, title = " ")+
geom_hline(yintercept=0, linetype="dashed")+scale_x_log10()+
theme_bw()
B.oe<-ggplot(b, aes(D50,Z.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
labs(y="Z-Beta", x="Dispersal (d50)", size=0.25, title=" ")+
theme_bw()+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")
S.cv.oe<-ggplot(b, aes(D50,CV.oe.S)) + geom_point(col="black", fill="grey55", shape=21,size=3)+
ylab("CV alpha real- CV alpha null")+
xlab("Dispersal (d50)")+
labs(title=Pondscapes_names[p])+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
B.cv.oe<-ggplot(b, aes(D50,CV.oe.Bet)) + geom_point(col="black", fill="grey55", shape=22,size=3)+
ylab("CV beta real- CV beta null")+
xlab("Dispersal (d50)")+
labs(title="")+scale_x_log10()+
geom_hline(yintercept=0, linetype="dashed")+
theme_bw()
supp_CV_plot[[p]]<- grid.arrange(S.cv.oe,B.cv.oe, ncol=2)
output_Figure3_data <- output_Figure3_data%>%bind_rows(
data.frame("Site" = "Rocha",
"Max_RatioAlpha"=max(f.hill.delta(d)),
"Dist_Max_RatioAlpha"=unique(d[which(f.hill.delta(d)==max(f.hill.delta(d)))]),
"Min_RatioBeta"=min(B.f.hill.delta(d)),
"Dist_Min_RatioBeta"=unique(d[which(B.f.hill.delta(d)==min(B.f.hill.delta(d)))]),
"Max_Z_Alpha"=max(b$Z.S),
"Dist_Max_Z_Alpha"=b$D50[which(b$Z.S==max(b$Z.S))],
"Min_Z_Alpha"=min(b$Z.S),
"Dist_Min_Z_Alpha"=b$D50[which(b$Z.S==min(b$Z.S))],
"Max_Z_Beta"=max(b$Z.Bet),
"Dist_Max_Z_Beta"=b$D50[which(b$Z.Bet==max(b$Z.Bet))],
"Min_Z_Beta"=min(b$Z.Bet),
"Dist_Min_Z_Beta"=b$D50[which(b$Z.Bet==min(b$Z.Bet))]
))
# Final plot assembly
final_plot_v1[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.oe, x = 0.85,y = 0, width = 0.15, height = 1)
final_plot_v2[[p]] <- ggdraw() +
draw_plot(S.o.Albera,x = 0, y = 0, width = 0.3, height = 1)+
draw_plot(B.o.Albera,x = 0.3, y = 0, width = 0.3, height = 1)+
draw_plot(legend, x = 0.6, y = 0.3, width = 0.1, height = 0.5)+
draw_plot(S.cv.oe, x = 0.7, y = 0, width = 0.15, height = 1)+
draw_plot(B.cv.oe, x = 0.85,y = 0, width = 0.15, height = 1)Figure 3 and Supplementary figures C and D printing
# Image printing
png(filename = "C:/Users/David CM/Dropbox/H2020/RETROMED/Frontiers_Diciembre_2022/Frontiers_Respuestas/Figures_Frontiers/Fig3_PondRetro_Deg_LGTBI.png",
width =1462*3, height =(length(Pondscapes_outputs)*(320*3)), units = "px",res = 300)
grid.arrange(final_plot_v1[[1]],
final_plot_v1[[2]],
final_plot_v1[[3]],
final_plot_v1[[4]],
final_plot_v1[[5]],
final_plot_v1[[6]],
nrow=6,ncol=1)
dev.off()
png(filename = "C:/Users/David CM/Dropbox/H2020/RETROMED/Frontiers_Diciembre_2022/Frontiers_Respuestas/Figures_Frontiers/Supplementary_Fig3_CV.png",
width =937*2, height =(length(Pondscapes_outputs)*(398*2)), units = "px",res = 300)
grid.arrange(supp_CV_plot[[1]],
supp_CV_plot[[2]],
supp_CV_plot[[3]],
supp_CV_plot[[4]],
supp_CV_plot[[5]],
supp_CV_plot[[6]],
nrow=6,ncol=1)
dev.off()
